library(conflicted) # 関数名の衝突を警告
conflicts_prefer( # 優先的に使う関数を指定
dplyr::filter(),
dplyr::select(),
dplyr::lag(),
)
library(tidyverse)
library(GGally)
library(broom) # 解析結果を tibble 形式に集約
library(gt) # 表の作成
library(gtsummary) # 分析結果の表の作成
library(ggfortify) # biplot表示のため第7講 主成分分析
評価と視覚化
準備
以下で利用する共通パッケージを読み込む.
寄与率・累積寄与率
分析結果の評価を行う関数としては以下がある.
base::summary(): 主成分負荷量や寄与率を表示broom::tidy(): 上記の内容を tibble 形式で取得ggfortify::autoplot(): 主成分得点による散布図
#' データフレーム 'toy_data' を分析
toy_pca <- toy_data |> select(計算対象の列) |> prcomp(必要なら標準化)
#' 主成分負荷量と寄与率を確認する
summary(toy_pca)
tidy(toy_pca, "d") # "u","v","d" で取得する情報を指定
#' 第1,2主成分得点で散布図を描く
autoplot(toy_pca, scale = 0) 以下では2つのデータセットを使用する
japan_social.csv の概要(再掲)
総務省統計局より取得した都道府県別の社会生活統計指標の一部
MASS::UScereal の概要
Nutritional and Marketing Information on US Cereals
The UScereal data frame has 65 rows and 11 columns. The data come from the 1993 ASA Statistical Graphics Exposition, and are taken from the mandatory F&DA food label. The data have been normalized here to a portion of one American cup.
詳細は ?MASS::UScereal を参照
標準化を行わずに分析する.
js_data <- read_csv("data/japan_social.csv") |>
mutate(Area = as_factor(Area))
js_data |> # 散布図.いくつかの変数は相関強いことがわかる
ggpairs(columns = 2:6, # 都道府県名・地方区分は除く
legend = c(2,1), # 2行1列のグラフから凡例を作成
lower = list(continuous = wrap("points", alpha = .5),
mapping = aes(colour = Area)))変数のばらつきに大きな違いがある.
js_data |> # 箱ひげ図.変数のばらつきに大きな違いがある
pivot_longer(where(is.double)) |> # 実数値をまとめる
mutate(name = as_factor(name)) |> # 列の順番どおりboxplotを並べる
ggplot(aes(x = name, y = value)) + # 既定値の name と value を利用
geom_boxplot(aes(fill = name), show.legend = FALSE) # 変数ごとに色を変える各変数の標本平均を0,不偏分散を1に規格化する.
js_data |> # 散布図.スケールは変わるが本質的には同じ図になる
mutate(across(where(is.double), \(x)c(scale(x)))) |>
ggpairs(columns = 2:6, # 都道府県名・地方区分は除く
legend = c(2,1), # 2行1列のグラフから凡例を作成
lower = list(continuous = wrap("points", alpha = .5),
mapping = aes(colour = Area)))標準化を行い,変数のばらつきを揃える.
js_data |> # 箱ひげ図.変数のばらつきをそろえる
mutate(across(where(is.double), \(x)c(scale(x)))) |>
pivot_longer(where(is.double)) |> # 実数値をまとめる
mutate(name = as_factor(name)) |> # 列の順番どおりboxplotを並べる
ggplot(aes(x = name, y = value)) + # 既定値の name と value を利用
geom_boxplot(aes(fill = name), show.legend = FALSE) # 変数ごとに色を変える標準化したデータを有効数字3桁で表示するには,例えば以下のようにすればよい.
js_data |> #
mutate(across(where(is.double), \(x)signif(c(scale(x)), digits = 3)))# A tibble: 47 × 7
Pref Forest Agri Ratio Land Goods Area
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <fct>
1 Hokkaido 0.425 4.63 0.979 -1.4 0.421 Hokkaido
2 Aomori 0.151 0.489 -0.512 -0.446 -0.274 Tohoku
3 Iwate 0.892 -0.159 -0.521 -0.776 -0.299 Tohoku
4 Miyagi -0.376 -0.361 -0.134 -1.1 0.993 Tohoku
5 Akita 0.599 -0.544 -0.614 -1.38 -0.48 Tohoku
6 Yamagata 0.479 0.205 -0.581 -0.574 -0.451 Tohoku
7 Fukushima 0.425 -0.734 -0.288 -1.08 -0.264 Tohoku
8 Ibaraki -2.04 0.691 0.0801 0.229 -0.123 Kanto
9 Tochigi -0.556 0.242 -0.269 -0.301 -0.127 Kanto
10 Gumma 0.151 0.994 -0.269 1.01 0.329 Kanto
# ℹ 37 more rows
問題
標準化の有無の違いで寄与率・累積寄与率がどのように異なるか確認しなさい.
prcomp(toy_data) # 標準化を行わない場合
prcomp(toy_data, scale. = TRUE) # 標準化を行う場合正式なオプション名は scale. であるが,sc = TRUE など他のオプションと区別できれば短縮表記も可能.
japan_social.csvの読み込み方の例js_data <- read_csv("data/japan_social.csv") |> mutate(Area = as_factor(Area))MASS::UScerealの整理の例glimpse(MASS::UScereal) # 各変数の属性を確認する. uc_data <- MASS::UScereal |> rownames_to_column(var = "product") |> # 行名の製品名を列に加える as_tibble()
base R の data.frame 型なので tibble 型に変換しておく.
解答例
MASS::UScereal の場合
元のデータは data.frame 形式なので,tibble に変換して整理しておく. また,適当な方法で視覚化をすることを推奨する.
uc_data <- MASS::UScereal |>
rownames_to_column(var = "product") |> as_tibble() 標準化なしの分析は以下のようになる.
uc_pca_raw <- uc_data |>
select(where(is.double)) |>
prcomp()
summary(uc_pca_raw)Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7
Standard deviation 203.19 98.6927 50.36206 6.82381 2.3718 1.48975 0.94740
Proportion of Variance 0.77 0.1817 0.04731 0.00087 0.0001 0.00004 0.00002
Cumulative Proportion 0.77 0.9517 0.99896 0.99983 0.9999 0.99998 0.99999
PC8
Standard deviation 0.55042
Proportion of Variance 0.00001
Cumulative Proportion 1.00000
uc_pca_raw$rotation PC1 PC2 PC3 PC4 PC5
calories -0.179207473 -0.1462324241 -0.96480715 0.045234798 -0.092696691
protein -0.011145612 0.0025160642 -0.01538123 -0.072061722 -0.169798607
fat -0.002912531 -0.0006278695 -0.01548280 0.077006299 -0.396417270
sodium -0.493173560 -0.8416992943 0.21967258 0.007757990 -0.001777641
fibre -0.027380174 0.0202681877 0.01023689 -0.075965904 -0.187070010
carbo -0.015056930 -0.0272462011 -0.10970338 -0.693413074 0.620733329
sugars -0.008609442 -0.0013826562 -0.04640333 0.707209145 0.620447068
potassium -0.850577043 0.5186488067 0.07824275 -0.005786195 0.012896751
PC6 PC7 PC8
calories -0.012034965 -0.0040178657 0.0696281364
protein -0.069385729 0.9054169859 -0.3755184499
fat 0.370507919 -0.3543044317 -0.7575401217
sodium -0.001439163 -0.0001605839 -0.0008059101
fibre -0.918767097 -0.2184131271 -0.2571603627
carbo 0.042042740 -0.0770477377 -0.3363935184
sugars -0.103827276 0.0320770971 -0.3175882779
potassium 0.032892502 -0.0016414154 0.0107594252
uc_pca_raw |>
tidy("d") |>
ggplot(aes(x = PC, y = percent)) + geom_bar(stat = "identity")uc_pca_raw |>
tidy("d") |>
ggplot(aes(x = PC, y = cumulative)) + geom_bar(stat = "identity")autoplot(uc_pca_raw, scale = 0,
data = uc_data, label = TRUE, label.label = "product") 標準化ありの分析は以下のようになる.
uc_pca <- uc_data |>
select(where(is.double)) |>
prcomp(scale. = TRUE)
summary(uc_pca)Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7
Standard deviation 2.0595 1.1595 1.0847 0.77916 0.70038 0.32184 0.17080
Proportion of Variance 0.5302 0.1681 0.1471 0.07589 0.06132 0.01295 0.00365
Cumulative Proportion 0.5302 0.6982 0.8453 0.92120 0.98252 0.99547 0.99911
PC8
Standard deviation 0.08427
Proportion of Variance 0.00089
Cumulative Proportion 1.00000
uc_pca$rotation PC1 PC2 PC3 PC4 PC5 PC6
calories -0.4095924 -0.40065656 0.1958834 0.007486492 -0.2032066 0.05255351
protein -0.4476405 0.17567026 0.0722067 -0.194640140 -0.1277506 -0.82959021
fat -0.2683607 -0.44417203 -0.2702646 -0.602181956 0.4973648 0.12339688
sodium -0.3474018 0.06764946 0.1195216 0.616762389 0.6907158 -0.03693360
fibre -0.3829156 0.46677617 -0.1976082 -0.097098382 -0.1367666 0.35604740
carbo -0.2882507 -0.20478232 0.6861866 -0.003106577 -0.2260198 0.29493081
sugars -0.1938910 -0.46258366 -0.5521416 0.455559106 -0.3701522 -0.03653264
potassium -0.4145359 0.36462256 -0.2330796 -0.045479335 -0.1054703 0.27809848
PC7 PC8
calories -0.35011143 0.68311726
protein 0.03187206 -0.14178644
fat 0.05108036 -0.17268441
sodium -0.03505495 -0.01931262
fibre -0.60446872 -0.27597329
carbo 0.23510504 -0.45909745
sugars 0.07382281 -0.30369583
potassium 0.66817773 0.32232225
uc_pca |>
tidy("d") |>
ggplot(aes(x = PC, y = percent)) + geom_bar(stat = "identity")uc_pca |>
tidy("d") |>
ggplot(aes(x = PC, y = cumulative)) + geom_bar(stat = "identity")autoplot(uc_pca_raw, scale = 0,
data = uc_data, label = TRUE, label.label = "product") 主成分分析の視覚化
主成分分析の分析結果を可視化する biplot 表示は, 関数 ggfortify::autoplot() にオプションを指定することで実行できる.
#' データフレームを分析する
toy_pca <-
toy_data |> select(計算対象の列を指定) |>
prcomp(必要なら標準化)
#' 第1と第2主成分を利用した biplot 表示
autoplot(toy_pca,
label = TRUE, # 各データのラベルを表示
loadings = TRUE, # 主成分負荷による成分方向の表示
loadings.label = TRUE) # 成分名の表示
#' 第2と第3主成分を利用した散布図
autoplot(toy_pca,
loadings = TRUE,
x = 2, y = 3) # 主成分の指定
#' パラメタ s を変更 (既定値は1)
autoplot(toy_pca,
scale = 0, # パラメタ s の指定
loadings = TRUE)base R には関数 biplot() が用意されている.
#' 第1と第2主成分を利用した散布図
biplot(toy_pca)
#' 第2と第3主成分を利用した散布図
biplot(toy_pca, choices = c(2,3))
#' パラメタ s を変更 (既定値は1)
biplot(toy_pca, scale=0)問題
それぞれのデータの主成分分析の結果を利用してバイプロットによる可視化を行いなさい.
- 標準化したデータでの主成分分析を行いなさい.
- 第1主成分と第2主成分でのバイプロットを描きなさい.
- 第2主成分と第3主成分でのバイプロットを描きなさい.
解答例
MASS::UScereal の場合
autoplot(uc_pca,
data = uc_data, colour = "mfr", # メーカー毎に色付け
shape = 19, size = 1,
label = TRUE, label.label = "product", label.repel = TRUE,
loadings = TRUE, loadings.colour = "orange",
loadings.label = TRUE, loadings.label.repel = TRUE,
loadings.label.colour = "brown", loadings.label.size = 4) +
theme(legend.position = "none")autoplot(uc_pca, x = 2, y = 3, # 第2 vs 第3
data = uc_data, colour = "mfr",
shape = 19, size = 1,
label = TRUE, label.label = "product", label.repel = TRUE,
loadings = TRUE, loadings.colour = "orange",
loadings.label = TRUE, loadings.label.repel = TRUE,
loadings.label.colour = "brown", loadings.label.size = 4) +
theme(legend.position = "none")第1,2主成分得点で散布図を描き,上と比較してみる.
augment(uc_pca, data = uc_data) |>
ggplot(aes(x = .fittedPC1, y = .fittedPC2, label = product, colour = mfr)) +
geom_point(shape = 19, size = 1) +
ggrepel::geom_text_repel(size = 3, max.overlaps = 40) +
theme(legend.position = "none")配置は変わらないが,座標軸が異なることがわかる. 関数 autoplot() において scale=0 とするとデータの座標は主成分得点となる.
options(ggrepel.max.overlaps = 40)
autoplot(uc_pca, scale = 0,
data = uc_data, colour = "mfr",
shape = 19, size = 1,
label = TRUE, label.label = "product", label.repel = TRUE, label.size = 3,
loadings = TRUE, loadings.colour = "orange",
loadings.label = TRUE, loadings.label.repel = TRUE,
loadings.label.colour = "brown", loadings.label.size = 4) +
theme(legend.position = "none")